Parameters

# Determining if doublet detecting and filtering will occur
if (isTRUE(params$doublet_filter)){
  dbl <- params$doublet_detect
  dbl_dtct=TRUE
} else {
  if (isFALSE(params$doublet_detect)) {
    dbl_dtct=FALSE
  } else {dbl_dtct=TRUE}
}

if (isTRUE(params$decontamX_filter)){
  dcx_dtct=TRUE
} else {
  if (isFALSE(params$decontamX_detect)) {
    dcx_dtct=FALSE
  } else {dcx_dtct=TRUE}
}


if (params$var_regress=="none"){
  var_regress=NULL
} else {
  var_regress <- params$var_regress
}

if ("cc" %in% params$var_regress) {
  cc_regress=TRUE
  var_regress=NULL
} else {
  cc_regress=FALSE
}

# This feature not yet supported.
if (!isFALSE(params$exp_metadata)) {
  cov_analysis <- TRUE
} else {
  cov_analysis <- FALSE
 }

if (params$parallel=="multicore") {
  bpparam <- MulticoreParam()
} else  {
  bpparam <- SerialParam()
}

Output Files

• The final seurat file will be saved to: ./results/QB_SingleCell_qQC/seurat_05Apr2022_0418UTC_10k.rds
  

seurat_raw <- readRDS(params$seurat_raw)

theme_sara <- function(){
  theme_bw(base_size=10)+
    theme(axis.text=element_text(color="black"),
          panel.background=element_rect(color="white"),
          strip.text = element_text(size=12),
          strip.background = element_rect(fill="white"))
}

theme_sara_90 <- function() {
  theme_bw(base_size=18)+
    theme(axis.text.x=element_text(angle=90,hjust = 1,vjust = 0.5),
          axis.text=element_text(color="black"),
          panel.background=element_rect(color="black"),
          strip.text = element_text(size=12),
          strip.background = element_rect(fill="white"))
}

sid <- params$sample_id
seurat_raw@meta.data$SampleID <- seurat_raw@meta.data[[sid]]

s.cols <- c(rcartocolor::carto_pal(n=4, "Teal"), 
            rcartocolor::carto_pal(n=4, "Emrld"), 
            rcartocolor::carto_pal(n=4, "BurgYl"), 
            rcartocolor::carto_pal(n=4, "PinkYl"), 
            rcartocolor::carto_pal(n=4, "Purp"))

num_samples <- length(levels(as.factor(seurat_raw@meta.data$SampleID)))
if (num_samples <= 20) {
  sample.cols <- sample(s.cols, size=num_samples, replace=FALSE)
} else {
  sample.cols <- colorRampPalette(s.cols)(num_samples)
}

names(sample.cols) <- levels(as.factor(seurat_raw@meta.data$SampleID))


heatmap.colfun <- microViz::heat_palette(
                  palette = "Lisbon",
                  breaks = 100,
                  range = c(0,1),
                  sym = FALSE,
                  rev = FALSE)

heatmap.colfun <- microViz::heat_palette(
                  palette = "Lisbon",
                  breaks = 100,
                  range = c(0,1),
                  sym = FALSE,
                  rev = FALSE)

dbl.cols <- c(`eSinglet`="snow2", `eDoublet`="#EE1289")
dcx.cols <- c(`Native`="azure2", `Ambient`="chartreuse1")
unikn_blue_red.cols <- c(rev(pal_karpfenblau[[5]]), "white", pal_bordeaux[[5]])


ref.cols <- list()

ref.cols[["SingleR_ImmGen"]]  <- c(`B cells` = "#F0A0FF", 
                                   `Basophils` = "#9370DB",
                                    `Eosinophils`= "#B452CD",
                                   `DC` = "#9DCC00", 
                                   `Epithelial cells` = "#87CEFA",  
                                   `Endothelial cells` = "#6495ED",
                                   `Fibroblasts` = "navy",
                                   `ILC` = "#191919", 
                                   `Macrophages` = "#1C8356", 
                                   `Mast cells` = "#E9Debb", 
                                   `Monocytes` = "#81C57A", 
                                   `Microglia` = "#EEB4B4",
                                   `Neutrophils` = "#b10da1",
                                   `NK cells` = "#8F7C00", 
                                   `NKT` = "#FEAF16", 
                                   `Others` = "gray40", 
                                   `Stem cells` = "red", 
                                   `Stromal cells` = "#F08080",
                                   `T cells` = "#FE902EFF", 
                                   `Tgd` = "#FFEE33")

ref.cols[["SingleR_HPCA"]] <- c(`Astocyte`="grey",
                                `B_cell`="#F0A0FF",
                                `BM`="chocolate4",
                                `BM & Prog.`="bisque3",
                                `CMP`="#DC8F71",
                                 `DC`="#9DCC00",  
                                `Endothelial_cells`="#6495ED",
                                `Epithelial_cells`="#87CEFA",
                                `Erythroblast`="red4",
                                `Fibroblasts`="navy",
                                `Gametocytes`="mediumpurple1",
                                `HSC`="red",
                                `Macrophage`="#1C8356",
                                `Monocyte`="#81C57A",
                                `Myelocyte`="green", 
                                `Neutrophils`="#b10da1",
                                `NK_cell`="#8F7C00",
                                `Osteoblasts`="cyan3",
                                `Platelets`="#A5DAC8",
                                `T_cells`="#FE902EFF",
                                `Others` = "gray40")
ref.cols[["SingleR_hpca"]] <- ref.cols[["SingleR_HPCA"]]                                

Distributions

Look for samples with skewed, abnormal or bimodal distributions

p1 <- RidgePlot(seurat_raw, group.by="SampleID", 
                features="nFeature_RNA", 
                same.y.lims=TRUE, sort=TRUE, log=TRUE)
p1 <- p1 + scale_fill_manual(values=sample.cols)
p1 <- p1 + scale_x_continuous(labels = scales::scientific)
p1 <- p1 + theme(axis.text=element_text(size=11))

p2 <- RidgePlot(seurat_raw, group.by="SampleID", 
                features="nCount_RNA", 
                same.y.lims=TRUE, sort=TRUE, log=TRUE)
p2 <- p2 + scale_fill_manual(values=sample.cols)
p2 <- p2 + scale_x_continuous(labels = scales::scientific, n.breaks=4)
p2 <- p2 + theme(axis.text=element_text(size=11))

mt.ul <- max(seurat_raw@meta.data$percent.mt)+5
p3 <- RidgePlot(seurat_raw, group.by="SampleID",
                features="percent.mt", y.max=50, sort=TRUE)
p3 <- p3 + scale_fill_manual(values=sample.cols)
p3 <- p3 + xlim(0, mt.ul)
p3 <- p3 + theme(axis.text=element_text(size=11))

cowplot::plot_grid(p1 + theme(legend.position="none", axis.title.y=element_blank()),
                           p2 + theme(legend.position="none", axis.title.y=element_blank()),
                           p3 + theme(legend.position="none", axis.title.y=element_blank()),
                           align = 'vh',
                           labels = c("A", "B", "C"),
                           hjust = -1,
                           nrow = 1)

Modality test

features <- c("nFeature_RNA", "nCount_RNA", "percent.mt")
names(features) <- features
modetest.list <- BiocParallel::bplapply(features, 
                   function(f) multimode::modetest(seurat_raw_20k@meta.data[[f]]),
                   BPPARAM=MulticoreParam())

pvalues_logical.list <- lapply(modetest.list, function(mt) mt[["p.value"]] <= 0.05)
                   
lapply(modetest.list, function(t) t)

## $nFeature_RNA
## 
##  Ameijeiras-Alonso et al. (2019) excess mass test
## 
## data:  seurat_raw_20k@meta.data[[f]]
## Excess mass = 0.034961, p-value < 2.2e-16
## alternative hypothesis: true number of modes is greater than 1
## 
## 
## $nCount_RNA
## 
##  Ameijeiras-Alonso et al. (2019) excess mass test
## 
## data:  seurat_raw_20k@meta.data[[f]]
## Excess mass = 0.0086293, p-value < 2.2e-16
## alternative hypothesis: true number of modes is greater than 1
## 
## 
## $percent.mt
## 
##  Ameijeiras-Alonso et al. (2019) excess mass test
## 
## data:  seurat_raw_20k@meta.data[[f]]
## Excess mass = 0.002569, p-value = 0.382
## alternative hypothesis: true number of modes is greater than 1

features.test <- features[unlist(pvalues_logical.list)]

modes.list <- BiocParallel::bplapply(features.test, 
                     function(f) multimode::locmodes(seurat_raw@meta.data[[f]],
                                                     mod0=2,
                                                     display=FALSE),
                     BPPARAM=bpparam)

antimode.list <- lapply(modes.list, 
                        function(m) round(m$locations[[2]], digits=2))

Inspect the distributions - if the antimode falls in the valley, it is potentially a good cutoff.

p1 <- ggplot(data=seurat_raw@meta.data, aes(x=nFeature_RNA)) 
p1 <- p1 + geom_density(fill="slateblue4") + labs(x="Number of genes")
if (isTRUE(pvalues_logical.list[["nFeature_RNA"]])) {
  p1 <- p1 + geom_vline(xintercept=antimode.list[["nFeature_RNA"]])
  p1 <- p1 + labs(caption=paste("antimode =", antimode.list[["nFeature_RNA"]]))
} else {
  p1 <- p1 + labs(caption="No antimode to display")
}
p1 <- p1 + theme(plot.caption=element_text(face="italic"))

p2 <- ggplot(data=seurat_raw@meta.data, aes(x=nCount_RNA))
p2 <- p2 + geom_density(fill="magenta4") + labs(x="Number of UMIs") + theme(axis.title.y=element_blank())
if (isTRUE(pvalues_logical.list[["nCount_RNA"]])) {
  p2 <- p2 + geom_vline(xintercept=antimode.list[["nCount_RNA"]])
  p2 <- p2 + labs(caption=paste("antimode =", antimode.list[["nCount_RNA"]]))
} else {
  p2 <- p2 + labs(caption="No antimode to display")
}
p2 <- p2 + theme(plot.caption=element_text(face="italic"))

mt.ul <- max(seurat_raw@meta.data$percent.mt)+2
p3 <- ggplot(data=seurat_raw@meta.data, aes(x=percent.mt))
p3 <- p3 + geom_density(fill="olivedrab3") + labs(x="% mitochondrial genes") + theme(axis.title.y=element_blank())
if (isTRUE(pvalues_logical.list[["percent.mt"]])) {
  p3 <- p3 + geom_vline(xintercept=antimode.list[["percent.mt"]])
  p3 <- p3 + labs(caption=paste("antimode =", antimode.list[["percent.mt"]]))
} else {
  p3 <- p3 + labs(caption="No antimode to display")
}
p3 <- p3 + theme(plot.caption=element_text(face="italic"))
p3 <- p3 +xlim(0, mt.ul)

cowplot::plot_grid(p1,
                   p2,
                   p3,
                   align = 'vh',
                   labels = c("A", "B", "C"),
                   hjust = -1,
                   nrow = 1)

Stats

Calculate median, median absolute deviation (mad), and outliers. Look for samples that are outliers.

median.list <- lapply(features, 
                       function(f) median(seurat_raw@meta.data[[f]], 
                                          na.rm=TRUE))

mad.list <- lapply(features,
                   function(f) mad(seurat_raw@meta.data[[f]],
                                      na.rm=TRUE))

outliers.list <- lapply(features,
                        function(f) scater::isOutlier(seurat_raw@meta.data[[f]], 
                                                      type="both"))
outliers.list <- lapply(outliers.list, function(o) attr(o, "thresholds"))

sample.median.df <- seurat_raw@meta.data %>%
                       group_by(SampleID) %>%
                       summarize(nFeature_RNA_median=median(nFeature_RNA, na.rm=TRUE),
                                 nCount_RNA_median=median(nCount_RNA, na.rm=TRUE),
                                 percent.mt_median=median(percent.mt, na.rm=TRUE))  

sample.median.df$nFeature_RNA_outlr <- ifelse(sample.median.df$nFeature_RNA_median <= median.list[["nFeature_RNA"]] - 3*mad.list[["nFeature_RNA"]] || sample.median.df$nFeature_RNA_median >= median.list[["nFeature_RNA"]] + 3*mad.list[["nFeature_RNA"]],
                                               "FLAG",
                                               "OK")



sample.median.df$nCount_RNA_outlr <- ifelse(sample.median.df$nCount_RNA_median <= median.list[["nCount_RNA"]] - 3*mad.list[["nCount_RNA"]] || sample.median.df$nFeature_RNA_median >= median.list[["nCount_RNA"]] + 3*mad.list[["nCount_RNA"]],
                                               "FLAG",
                                               "OK")

sample.median.df$percent.mt_outlr <- ifelse(sample.median.df$percent.mt_median >= median.list[["percent.mt"]] + 3*mad.list[["percent.mt"]],
                                               "FLAG",
                                               "OK")

kable(sample.median.df, booktabs = TRUE) %>%
  kable_styling(font_size = 12)

All metrics

keep_remove.cols <- c(`Keep`="lavenderblush3",
                      `Remove`="#DC143C")

p1 <- ggplot(data=seurat_raw@meta.data, aes(x=nCount_RNA, y=nFeature_RNA, color=Filter))
p1 <- p1 + geom_point() + geom_hline(yintercept=nFeature_RNA.ll) + geom_hline(yintercept=nFeature_RNA.ul) + geom_vline(xintercept=nCount_RNA.ul)
p1 <- p1 + theme_sara()
p1 <- p1 + scale_color_manual(values=keep_remove.cols)

p2 <- ggplot(data=seurat_raw@meta.data, aes(x=percent.mt, y=nFeature_RNA, color=Filter))
p2 <- p2 + geom_point() + geom_hline(yintercept=nFeature_RNA.ll) + geom_hline(yintercept=nFeature_RNA.ul) + geom_vline(xintercept = percent.mt.ul)
p2 <- p2 + theme_sara()
p2 <- p2 + scale_color_manual(values=keep_remove.cols)

p3 <- ggplot(data=seurat_raw@meta.data, aes(x=percent.mt, y=nCount_RNA, color=Filter))
p3 <- p3 + geom_point() + geom_vline(xintercept = percent.mt.ul) + geom_hline(yintercept=nCount_RNA.ul)
p3 <- p3 + geom_hline(yintercept=nCount_RNA.ll)
p3 <- p3 + theme_sara()
p3 <- p3 + scale_color_manual(values=keep_remove.cols)

legend <- cowplot::get_legend(p1 + theme(legend.box.margin = margin(0, 0, 0, 12)))

prow <- cowplot::plot_grid(p1 + theme(legend.position="none"),
                           p2 + theme(legend.position="none"),
                           p3 + theme(legend.position="none"),
                           align = 'vh',
                           labels = c("A", "B", "C"),
                           hjust = -1,
                           nrow = 1)
cowplot::plot_grid(prow, legend, rel_widths = c(4, .6))

Idents(seurat_raw) <- seurat_raw@meta.data$orig.ident
seurat <- subset(seurat_raw, subset=Filter=="Keep")

Inspect samples after filtering

Look for outlier samples to make note of.

meta.df <- as.data.table(seurat@meta.data)
summary.df <- meta.df %>% 
              dplyr::group_by(SampleID) %>% 
              dplyr::summarise(n=n(),
                               median_nFeature=median(nFeature_RNA, na.rm=TRUE),
                               medan_nCount=median(nCount_RNA, na.rm=TRUE),
                               median_percentMT=median(percent.mt, na.rm=TRUE))
setDT(summary.df)
summary.melt.df <- data.table::melt(summary.df, id.vars="SampleID")
setorder(summary.df, n)
x.ord <- summary.df$SampleID

p1 <- ggplot(data=summary.df, aes(x=SampleID, y=n, fill=SampleID))
p1 <- p1 + geom_bar(stat="identity", color="black")
p1 <- p1 + scale_fill_manual(values=sample.cols)
p1 <- p1 + theme_sara_90() + theme(axis.text=element_text(size=9), axis.title.x = element_blank(), axis.title.y=element_text(size=11))
p1 <- p1 + scale_x_discrete(limits=x.ord)
p1 <- p1 + labs(y="Number of cells")

p2 <- ggplot(data=meta.df, aes(x=SampleID, y=nFeature_RNA, fill=SampleID))
p2 <- p2 + geom_boxplot()
p2 <- p2 + scale_fill_manual(values=sample.cols)
p2 <- p2 + theme_sara_90() + theme(axis.text=element_text(size=9), axis.title.x = element_blank(), axis.title.y=element_text(size=11))
p2 <- p2 + scale_x_discrete(limits=x.ord)
p2 <- p2 + labs(y="Number of Genes (nFeature_RNA)")

p3 <- ggplot(data=meta.df, aes(x=SampleID, y=nCount_RNA, fill=SampleID))
p3 <- p3 + geom_boxplot()
p3 <- p3 + scale_fill_manual(values=sample.cols)
p3 <- p3 + theme_sara_90() + theme(axis.text=element_text(size=9), axis.title.x = element_blank(), axis.title.y=element_text(size=11))
p3 <- p3 + scale_x_discrete(limits=x.ord)
p3 <- p3 + labs(y="Number of UMIs (nCount_RNA)")

p4 <- ggplot(data=meta.df, aes(x=SampleID, y=percent.mt, fill=SampleID))
p4 <- p4 + geom_boxplot()
p4 <- p4 + scale_fill_manual(values=sample.cols)
p4 <- p4 + theme_sara_90() + theme(axis.text=element_text(size=9), 
                                   axis.title.x = element_blank(), axis.title.y=element_text(size=11))
p4 <- p4 + scale_x_discrete(limits=x.ord)
p4 <- p4 + labs(y="Percentage MT")

#legend <- cowplot::get_legend(p2 + theme(legend.box.margin = margin(0, 0, 0, 12)))

cowplot::plot_grid(p1 + theme(legend.position="none"),
                   p2 + theme(legend.position="none"),
                   p3 + theme(legend.position="none"),
                   p4 + theme(legend.position="none"),
                   align = 'vh',
                   labels = c("A", "B", "C", "D"),
                   hjust = -1,
                   nrow = 2)

Cell cycle regression

If cell cycle was previously identified as having an effect on the clustering and the differences are not biological, the difference between G2M and S phase scores will be calculated and used as the variable to regress. For most cases, it is not recommended to regress out cell cycle. Other variables will also be regressed out at this time.

seurat@meta.data$cc_dif <- seurat@meta.data$S.Score - seurat@meta.data$G2M.Score

var_regress <- gsub("cc", "cc_dif", params$var_regress)

seurat <- ScaleData(object = seurat,
                        vars.to.regress=var_regress,
                        verbose = TRUE)

Prepare reference data

The Human Primary Cell Atlas (HPCA, human) and ImmGen (mouse) are both good choices. Mouse data can also be annotated with HPCA. HPCA is written as HPCA when the data is from human or when human reference gene names are converted to mouse gene names via biomart. HPCA is denoted as hpca when the reference data is converted to mouse names by changing to sentence case.

ref.data.list <- list()

if (params$organism=="human") {
  if ("HPCA" %in% params$anno_db) {
    ref.data.list[["HPCA"]] <- celldex::HumanPrimaryCellAtlasData(ensembl=FALSE) 
  }
}

if (params$organism=="mouse") {
  if ("HPCA" %in% params$anno_db) {
      #human <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
      #mouse <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")
      ref.data <- celldex::HumanPrimaryCellAtlasData(ensembl=FALSE)  
      hpa_human.genes <- rownames(assay(ref.data))

      human_to_mouse.df <- readRDS(params$biomart_rds)
      #human_to_mouse.df <- getLDS(attributes=c("external_gene_name"),
      #                            filters="external_gene_name",
      #                            values = hpa_human.genes, 
      #                            mart = human,
      #                            attributesL = c("external_gene_name"), martL = mouse)
    
      assay.data <- as.data.frame(assay(ref.data))
      assay.data$Gene.name <- rownames(assay.data)
      assay.mdata <- merge(assay.data, human_to_mouse.df, by="Gene.name")
      assay.mdata$Gene.name <- NULL
      assay.mdata <- assay.mdata[!duplicated(assay.mdata$Gene.name.1),]
      assay.mdata <- tibble::remove_rownames(assay.mdata)
      assay.mdata <- tibble::column_to_rownames(assay.mdata, "Gene.name.1")

      # Lose about 3k marker genes out of 19k
      dim(assay.mdata)
      dim(assay.data)

      # over-write ref data
      h2m_ref.data <- SummarizedExperiment(assays=SimpleList(logcounts=as.matrix(assay.mdata)),
                                     rowData=rownames(assay.mdata),
                                     colData=colData(ref.data),
                                     metadata=metadata(ref.data))
      ref.data.list[["HPCA"]] <- h2m_ref.data
  }    
  if ("hpca" %in% params$anno_db) {
      ref.data <- celldex::HumanPrimaryCellAtlasData(ensembl=FALSE)  
      hpa_human.genes <- rownames(assay(ref.data))
      
      assay.data <- as.data.frame(assay(ref.data))
      assay.data$Gene.name <- rownames(assay.data)
      assay.data$Gene.name <- stringr::str_to_sentence(assay.data$Gene.name)
      rownames(assay.data) <- assay.data$Gene.name
      assay.data$Gene.name <- NULL
      
      # overwrite reference data
      h2m_ref.data <- SummarizedExperiment(assays=SimpleList(logcounts=as.matrix(assay.data)),
                                           rowData=rownames(assay.data),
                                           colData=colData(ref.data),
                                           metadata=metadata(ref.data))
      ref.data.list[["hpca"]] <- h2m_ref.data
  }
  if ("ImmGen" %in% params$anno_db) {
  ref.data.list[["ImmGen"]] <- celldex::ImmGenData(ensembl=FALSE)  
  }
}

SingleR annotation

The differential expression method of Wilcoxon is used for training the classifier, which often provides a better annotation over the classic method. See links below. Spearman’s correlation is used for classification.

The Single R Book

SingleR marker detection

cells.list <- BiocParallel::bplapply(ref.data.list,
                                     function(ref) SingleR(test=seurat@assays[["RNA"]]@data,
                                                           ref=ref, 
                                                           labels=ref$label.main, 
                                                           fine.tune = TRUE, 
                                                           de.method="wilcox", 
                                                           de.n=30, 
                                                           BPPARAM=bpparam)) 
### Add to seurat
for (ann in names(cells.list)) {
  ann_name <- paste0("SingleR_", ann)
  ann_name_pruned <- paste0(ann_name, "_pruned")
  seurat@meta.data[[ann_name]] <- cells.list[[ann]]$labels
  seurat@meta.data[[ann_name_pruned]] <- cells.list[[ann]]$pruned.labels
}

num_ann <- length(cells.list)
nrow <- round((num_ann+.1)/2, digits=0)
fh <- nrow*6

Assess annotation

rand_cells <- runif(500, min=1, max=500)
plot_shmap <- function(anno_db, cells) {
  p <- plotScoreHeatmap(cells, 
                        cells.use=rand_cells, 
                        silent=TRUE,
                        treeheight_row=0,
                        fontsize=9,
                        main=anno_db)
  return(p[[4]])
}

shmap.plots <- sapply(names(cells.list), 
                      function(a) plot_shmap(anno_db=a,
                                             cells=cells.list[[a]]),
                      simplify=FALSE, USE.NAMES=TRUE)



cowplot::plot_grid(plotlist=shmap.plots,
                   align="vh",
                   labels="AUTO",
                   label_size=12,
                   ncol=2,
                   nrow=nrow)

plot_delta <- function(anno_db, cells) {
  p <- plotDeltaDistribution(cells, 
                      ncol=5,
                      size=0.5) +
                      theme_sara()+ theme(strip.text = element_text(size=9),
                                          legend.position="bottom",
                                          axis.text.x=element_blank())+
                      labs(y="Delta score", x="", title=anno_db)
  return(p)
}

delta.plots <- sapply(names(cells.list), 
                      function(a) plot_delta(anno_db=a,
                                             cells=cells.list[[a]]),
                      simplify=FALSE, USE.NAMES=TRUE)

cowplot::plot_grid(plotlist=delta.plots,
                   align="vh",
                   labels="AUTO",
                   label_size=12,
                   ncol=2,
                   nrow=nrow)

pa <- params$pref_ann
pref_ann <- paste0("SingleR_", pa)

prune <- pruneScores(cells.list[[pa]],
                   nmads = 3,
                   min.diff.med = -Inf,
                   min.diff.next = 0,
                   get.thresholds = FALSE)

prune <- which(prune=="TRUE")
keep <- cells.list[[pa]][-prune,]
seurat <- seurat[,-prune]

SingleR_names <- paste0("SingleR_", names(cells.list))
names(SingleR_names) <- SingleR_names

collapse_other <- function(seurat, ref_label, freq.ll, cells.list) {
  seurat@meta.data[[ref_label]] <- factor(seurat@meta.data[[ref_label]])
  freq<-data.frame(janitor::tabyl(seurat@meta.data[[ref_label]], sort = TRUE))
  colnames(freq)<-c("cells","n","proportion")
  freq$percent<-freq$proportion*100
  freq<-freq[order(freq$n, decreasing = TRUE),]
  others<-as.vector(freq$cells[which(freq$n<freq.ll)])
  seurat@meta.data[[ref_label]] <- forcats::fct_collapse(seurat@meta.data[[ref_label]],
                                           Others=c(others))
  return(seurat@meta.data[[ref_label]])
}

for (a in SingleR_names) {
  seurat@meta.data[[a]] <- collapse_other(seurat,
                                                  ref_label=a,
                                                  freq.ll=params$ann_collapse_n,
                                                  cells.list=cells.list)
}

if ("HPCA" %in% params$anno_db) {
  seurat@meta.data$SingleR_HPCA <- as.factor(seurat@meta.data$SingleR_HPCA)
  seurat@meta.data$SingleR_HPCA <- forcats::fct_collapse(seurat@meta.data$SingleR_HPCA, 
                                                         B_cell=c("B_cell", "Pro-B_cell_CD34+", "Pre-B_cell_CD34-"),
                                                         Others=c("Others", "Embryonic_stem_cells", "GMP",
                                                                  "Hepatocytes", "iPS_cells", "Keratinocytes",
                                                                  "MSC", "Neuroepithelial_cell", "Neurons"),
                                                         HSC=c("HSC_-G-CSF", "HSC_CD34+"),
                                                         Fibroblasts=c("Fibroblasts", "Chondrocytes",
                                                                       "Smooth_muscle_cells",
                                                                       "Tissue_stem_cells"),
                                                         Myelocyte=c("Myelocyte", "Pro-Myelocyte"))
    
}
if ("hpca" %in% params$anno_db) {
 seurat@meta.data$SingleR_hpca <- as.factor(seurat@meta.data$SingleR_hpca)
  seurat@meta.data$SingleR_hpca <- forcats::fct_collapse(seurat@meta.data$SingleR_hpca, 
                                                         B_cell=c("B_cell", "Pro-B_cell_CD34+", "Pre-B_cell_CD34-"),
                                                         Others=c("Others", "Embryonic_stem_cells", "GMP",
                                                                  "Hepatocytes", "iPS_cells", "Keratinocytes",
                                                                  "MSC", "Neuroepithelial_cell", "Neurons"),
                                                         HSC=c("HSC_-G-CSF", "HSC_CD34+"),
                                                         Fibroblasts=c("Fibroblasts", "Chondrocytes",
                                                                       "Smooth_muscle_cells",
                                                                       "Tissue_stem_cells"),
                                                          Myelocyte=c("Myelocyte", "Pro-Myelocyte"))
}
if ("ImmGen" %in% params$anno_db) {
    seurat@meta.data$SingleR_ImmGen <- as.character(seurat@meta.data$SingleR_ImmGen)
    seurat@meta.data$SingleR_ImmGen <- if_else(seurat@meta.data$SingleR_ImmGen=="B cells, pro", 
                                               "B cells",
                                                seurat@meta.data$SingleR_ImmGen)
}

dat <- as.data.table(seurat@meta.data)
ann <- paste0("SingleR_", params$pref_ann)

dt <- dat[, .N, by=ann] 
dt <- dt[order(dt[[ann]]),]


if (ann %chin% names(ref.cols)) {
  pref_ann.cols <- ref.cols[[ann]]
} else {
  pref_ann.cols <- colorRampPalette(pals::watlington(n=16))(nrow(dt))
  names(pref_ann.cols) <- dt[[ann]]
}

p <- ggplot(dt, aes(x="", y=N, fill=.data[[ann]]))
p <- p + geom_bar(stat="identity", width=1, color="white")
p <- p + coord_polar("y", start=0)
p <- p + scale_fill_manual(values=pref_ann.cols)
p <- p + theme_void()  
p

Elbow plot

The elbow plot gives us an idea of how many PCs will be significant. When the plot plateaus, that is the point of diminishing returns.

dat <- data.frame(`Standard_deviation`=seurat@reductions[["pca"]]@stdev,
                  PC=c(1:length(seurat@reductions[["pca"]]@stdev)))
p <- ggplot(data=dat, aes(x=PC, y=Standard_deviation))
p <- p + geom_bar(stat="identity", fill="#6AAAB7", color=NA)
p <- p + theme_sara()
p <- p + labs(y="Standard deviation")
p

Dimension loadings

The dimension loadings give us an idea of the top genes contributing to each PC.

plot.list <- VizDimLoadings(seurat, dims=1:9, ncol=3, combine=FALSE, col=NA)
for (i in 1:length(plot.list)) {
  pc <- paste("PC", i, sep="_")
  #scale_lim <- max(a)
  p <- plot.list[[i]]
  xscale_lim <- max(abs(p$data[[pc]]))
  p <- p + geom_bar(stat="identity", aes(fill=.data[[pc]])) 
  p<- p + scale_fill_gradient2(low=unikn_blue_red.cols[1], 
                               mid="white",
                               high=unikn_blue_red.cols[3],
                               limits=c(-xscale_lim, xscale_lim))
  p <- p + theme_sara() 
  p <- p + theme(axis.title = element_blank(),
                 axis.text.x=element_text(angle=90))
  plot.list[[i]] <- p
}

cowplot::plot_grid(plotlist=plot.list,
                   align="vh",
                   labels="AUTO",
                   label_size=12,
                   ncol=3)

How well do clusters agree with annotation?

We calculate the adjusted rand index (ARI) for clusters vs annotations, which ranges from 0 to 1. Larger values indicate a better fit. Ideally, we would hope for an ARI above at least 0.5, although this can vary with the single-cell data and availability of suitable references.

cluster.names <- names(seurat@meta.data)[grep("Leiden|Louvain", names(seurat@meta.data))]
names(cluster.names) <- cluster.names
ARI.list <- lapply(cluster.names, 
                   function(c) lapply(SingleR_names, 
                                      function(a) mclust::adjustedRandIndex(seurat@meta.data[[a]],
                                                                            seurat@meta.data[[c]])))
ARI.list <- unlist(ARI.list, recursive=FALSE)
ARI.dt <- data.table(ID=names(ARI.list),
                     ARI=unlist(ARI.list))

ARI.dt <- ARI.dt[, c("Cluster") := tstrsplit(ID, "_", keep=1)]
ARI.dt <- ARI.dt[, c("Resolution", "Annotation") := tstrsplit(ID, "_", keep=c(2,3))]
ARI.dt$Resolution <- gsub(".SingleR", "", ARI.dt$Resolution)
ARI.dt$Annotation <- paste0("SingleR_", ARI.dt$Annotation)
ARI.dt$cluster_res <- paste(ARI.dt$Cluster, ARI.dt$Resolution, sep="_")
num_ann <- length(levels(as.factor(ARI.dt$Annotation)))
fw <- (num_ann*2.5)+1
nr <- round((num_ann+.1)/2, digits=0)
fh <- (nr*2.5)+.5

# karpfenblau2 "#B4BCD6", "#8290BB",  "#586BA4", "#3E5496","#324376
res.cols <- colorRampPalette(pal_karpfenblau)(length(params$cluster_resolution))
names(res.cols) <- levels(as.character(params$cluster_resolution))
p1 <- ggplot(data=ARI.dt, aes(x=Cluster, y=ARI, fill=as.character(Resolution)))
p1 <- p1 + geom_bar(stat="identity", color="white", position=position_dodge(preserve="total")) + scale_fill_manual(name="Resolution", values=res.cols)
p1 <- p1 + ylim(0,1)
p1 <- p1 + facet_wrap(~Annotation, scale="free_x", nrow=nr)
p1 <- p1 + theme_sara() + theme(strip.text = element_text(size=11))
p1

# Pick top ARI regardless of preferred annotation
setorder(ARI.dt, -ARI)
ARI.dt2 <- ARI.dt[, c("Annotation", "Cluster", "Resolution", "ARI"), with=FALSE]
ARI.dt2$ARI <- round(ARI.dt2$ARI, digits=4)
kable(ARI.dt2, booktabs = TRUE) %>%
  kable_styling(font_size = 12)

top_ann <- ARI.dt$Annotation[1]
top_clust_res <- ARI.dt$cluster_res[1]

if (top_ann != pref_ann) {
  pref.ARI.dt <- ARI.dt[Annotation==pref_ann]
  setorder(pref.ARI.dt, -ARI)
  pref_clust_res <- pref.ARI.dt$cluster_res[1]
  poor_choice=TRUE
} else {
  pref_clust_res <- top_clust_res
  poor_choice=FALSE
}

Comparison of UMAP by cluster vs annotation

UMAP clusters
Inspecting the UMAP This will give us an idea of how many clusters and how well separated. Look for the cells to stay within their cluster.

top_clust.cols <- colorRampPalette(carto_pal(n=12, "Vivid"))(length(levels(as.factor(seurat@meta.data[[top_clust_res]]))))
names(top_clust.cols) <- levels(as.factor(seurat@meta.data[[top_clust_res]]))


p <- DimPlot(seurat, group.by=top_clust_res) 
p <- p + scale_color_manual(values=top_clust.cols)
p <- p + theme_sara()
p

if (top_ann %chin% names(ref.cols)) {
  top_ann.cols <- ref.cols[[top_ann]]
} else {
   top_ann.cols <- colorRampPalette(pals::alphabet2(n=26))(nrow(dt))
   names(top_ann.cols) <- levels(as.factor(seurat@meta.data[[top_ann]]))

}

p <- DimPlot(seurat, group.by=top_ann) 
p <- p + scale_color_manual(values=top_ann.cols)
p <- p + theme_sara()
p

p <- DimPlot(seurat, group.by=pref_ann) 
p <- p + scale_color_manual(values=pref_ann.cols)
p <- p + theme_sara()
p <- p + theme(plot.title=element_blank())
p

dat <- seurat@meta.data
p <- ggplot(data=dat, aes(x=.data[[top_clust_res]], fill=.data[[top_ann]]))
p <- p + geom_bar(position="fill", color="black")
p <- p + theme_sara_90()
p <- p + scale_fill_manual("", values=top_ann.cols)
p <- p + labs(y="Proportion")
p <- p + theme(axis.title.x=element_blank(),
               axis.text=element_text(size=9),
               legend.text=element_text(size=9),
               legend.title=element_text(size=9),
               axis.title.y=element_text(size=10),
               strip.text=element_text(size=10),
               legend.position="bottom")
p

 dat <- seurat@meta.data
p <- ggplot(data=dat, aes(x=.data[[pref_clust_res]], fill=.data[[pref_ann]]))
p <- p + geom_bar(position="fill", color="black")
p <- p + theme_sara_90()
p <- p + scale_fill_manual("", values=pref_ann.cols)
p <- p + labs(y="Proportion")
p <- p + theme(axis.title.x=element_blank(),
               axis.text=element_text(size=9),
               legend.text=element_text(size=9),
               legend.title=element_text(size=9),
               axis.title.y=element_text(size=10),
               strip.text=element_text(size=10),
               legend.position="bottom")
p

Cluster markers

If the cluster and annotations don’t match up well, it can be helpful to look at the markers that differentiate the clusters. Cell cycle markers, for example, may be problematic.

Idents(seurat) <- seurat@meta.data[[top_clust_res]]
markers <- FindAllMarkers(seurat, 
                          only.pos = TRUE,
                          min.pct = 0.3, 
                          logfc.threshold = .5,
                          max.cells.per.ident=5000,
                          random.seed = 888)

markers %>%
    group_by(cluster) %>%
    top_n(n = 5, wt = avg_log2FC) -> top5
marker.genes <- top5$gene
seurat.hm <- subset(seurat, downsample=1000)
mat<- seurat.hm[["RNA"]]@data[marker.genes, ] %>% as.matrix()
mat<- log2(mat+1)
#mat <- seurat.hm[["RNA"]]@scale.data[marker.genes,] %>% as.matrix() # use this and dont log to get DoHeatmap style

# cell colors
cell_anno<- seurat.hm@meta.data[[pref_ann]]
cell.levels <- levels(droplevels(as.factor(cell_anno)))
c.cols <- pref_ann.cols[names(pref_ann.cols) %in% cell.levels]

# cluster colors
clust_anno <- as.character(seurat.hm@meta.data[[top_clust_res]])
clust_levels <- levels((droplevels(as.factor(clust_anno))))
clust.cols <- colorRampPalette(carto_pal(n=12, name="Vivid"))(length(clust_levels))
names(clust.cols) <- clust_levels
  

# cluster rowcolors
clust_anno2 <- as.character(top5$cluster)
clust_levels <- levels((droplevels(as.factor(clust_anno2))))
clust.cols2 <- clust.cols


ta <- ComplexHeatmap::HeatmapAnnotation(`Cluster`=clust_anno,
                                        `Cell type`=cell_anno,
                                         col=list(`Cell type`=c.cols,
                                                  `Cluster`=clust.cols))

ra <- ComplexHeatmap::rowAnnotation(`Cluster`=clust_anno2,
                                     col=list(`Cluster`=clust.cols2),
                                     show_legend=FALSE)

heatmap2.colfun <- microViz::heat_palette(palette = "Cividis",
                                          breaks = 100,
                                          range = c(0,max(mat)),
                  sym = FALSE,
                  rev = FALSE)

ComplexHeatmap::Heatmap(mat, name = "Expression",  
                        cluster_columns = FALSE,
                        cluster_rows=FALSE,
                        column_split=seurat.hm@meta.data[[top_clust_res]],
                        row_split=top5$cluster,
                        row_names_gp = grid::gpar(fontsize = 11),
                        column_title=character(0),
                        column_gap = unit(0.5, "mm"),
                        col = heatmap2.colfun,
                        top_annotation = ta,
                        show_column_names = FALSE,
                        left_annotation = ra)

Is there an effect of cell cycle?

There may be some structure in the UMAP due to cell cycle, it is also a good idea to check the PCs and cluster markers.

p <- DimPlot(seurat, group.by="Phase") 
p <- p + scale_color_manual(values=cell_cycle.cols)
p <- p + theme_sara()
p

tab<-prop.table(table(seurat@meta.data$Phase, seurat@meta.data[[top_clust_res]]),2)

nr <- length(levels(as.factor(seurat@meta.data$Phase)))
nc <- length(levels(as.factor(seurat@meta.data[[top_clust_res]])))


ComplexHeatmap::Heatmap(tab,
                        cluster_columns = TRUE,
                        show_column_dend = TRUE,
                        row_names_gp = grid::gpar(fontsize = 12),
                        column_names_gp = grid::gpar(fontsize=12),
                        cluster_rows = TRUE,
                        show_row_dend = FALSE,
                        col = heatmap.colfun,
                        height = unit(5, "mm")*nr,
                        width = unit(5, "mm")*nc,
                        name="Proportion")

tab<-prop.table(table(seurat@meta.data$Phase, seurat@meta.data[[top_ann]]),2)

nr <- length(levels(as.factor(seurat@meta.data$Phase)))
nc <- length(levels(as.factor(seurat@meta.data[[top_ann]])))


ComplexHeatmap::Heatmap(tab,
                        cluster_columns = TRUE,
                        show_column_dend = TRUE,
                        row_names_gp = grid::gpar(fontsize = 12),
                        column_names_gp = grid::gpar(fontsize=12),
                        cluster_rows = TRUE,
                        show_row_dend = FALSE,
                        col = heatmap.colfun,
                        height = unit(5, "mm")*nr,
                        width = unit(5, "mm")*nc,
                        name="Proportion")

Cell metrics effects

First we will look at the UMAP for visual inspection. Extremely low or high UMI/genes in a cluster or high percent.mt in a cluster may be problematic.

feature_plot <- function(seurat, feature) {
  p <- ggplot(data=as.data.frame(seurat@reductions$umap@cell.embeddings), 
              aes(x=UMAP_1, y=UMAP_2, 
                  color=seurat@meta.data[[feature]]))
  p <- p + geom_point(size=0.1)
  p <- p + theme_sara()
  if(feature == "percent.mt") {
      p <- p + scale_color_viridis_c(option="rocket", 
                                 direction=-1)
  } else {
    p <- p + scale_color_viridis_c(option="rocket", 
                                 direction=-1,
                                 labels=scales::label_scientific())
  }
  p <- p + theme(axis.title.y.right = element_blank(),
                 legend.position="right",
                 legend.title=element_blank())
  p <- p + labs(title=feature)
return(p)
}

p.list <- sapply(names(features), 
                 function(feature) feature_plot(seurat, feature),
                 simplify=FALSE, USE.NAMES=TRUE)

cowplot::plot_grid(plotlist=p.list,
                   align="vh",
                   labels="AUTO",
                   label_size=12,
                   ncol=2,
                   nrow=2)

p1 <- RidgePlot(seurat, group.by=top_clust_res, features="nFeature_RNA", same.y.lims=TRUE, sort=TRUE) 
p1 <- p1 + scale_fill_manual(values=top_clust.cols)
p1 <- p1 + theme(legend.position="none", 
                 axis.title.y=element_blank(),
                 axis.text.x=element_text(size=11))
p1 <- p1 + scale_x_continuous(labels = scales::scientific)

p2 <- RidgePlot(seurat, group.by=top_clust_res, features="nCount_RNA", same.y.lims=TRUE, sort=TRUE, log=TRUE)
p2 <- p2 + scale_fill_manual(values=top_clust.cols)
p2 <- p2 + theme(legend.position="none", 
                 axis.text=element_text(size=11),
                 axis.title.y=element_blank())
p2 <- p2 + scale_x_continuous(labels = scales::scientific)

mt.ul <- max(seurat@meta.data$percent.mt)+5
p3 <- RidgePlot(seurat, group.by=top_clust_res, features="percent.mt", y.max=50, sort=TRUE)
p3 <- p3 + scale_fill_manual(values=top_clust.cols)
p3 <- p3 + xlim(0, mt.ul)
p3 <- p3 + theme(legend.position="none", 
                 axis.title.y=element_blank(),
                 axis.text=element_text(size=11))

prow <- cowplot::plot_grid(p1,
                           p2,
                           p3,
                           align = 'vh',
                           labels = "AUTO",
                           hjust = -1,
                           nrow = 1)
prow

p1 <- RidgePlot(seurat, group.by=top_ann, features="nFeature_RNA", same.y.lims=TRUE, sort=TRUE) 
p1 <- p1 + scale_fill_manual(values=top_ann.cols)
p1 <- p1 + theme(legend.position="none", axis.title.y=element_blank())

p2 <- RidgePlot(seurat, group.by=top_ann, features="nCount_RNA", same.y.lims=TRUE, sort=TRUE, log=TRUE)
p2 <- p2 + scale_fill_manual(values=top_ann.cols)
p2 <- p2 + theme(legend.position="none", axis.title.y=element_blank())

p3 <- RidgePlot(seurat, group.by=top_ann, features="percent.mt", y.max=50, sort=TRUE)
p3 <- p3 + scale_fill_manual(values=top_ann.cols)
p3 <- p3 + theme(legend.position="none", axis.title.y=element_blank())
p3 <- p3 + xlim(0, percent.mt.ul+2)

prow <- cowplot::plot_grid(p1,
                           p2,
                           p3,
                           align = 'vh',
                           labels = "AUTO",
                           hjust = -1,
                           nrow = 1)
prow